home *** CD-ROM | disk | FTP | other *** search
/ AOL File Library: 2,801 to 2,900 / aol-file-protocol-4400-2801-to-2900.zip / AOLDLs / C++ Files Library / Graphic Gems I, II & III (C_C++) / Graphics Gems C Code.sea / GemsIII / accurate_scan / fixpoint.c < prev    next >
C/C++ Source or Header  |  1992-06-16  |  5KB  |  227 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "fixpoint.h"
  5.  
  6. /* A non-optimized fixpoint pkg, with some overflow checking. */
  7.  
  8. #define TwosComplement(_x_)  (-(_x_))
  9.  
  10. int fp_verbose;
  11. int fp_error = 0;        /* will contain error indicators for fp_multiply */
  12. int fp_print_error = 1;  /* if 1, print overflow errors */
  13.  
  14.  
  15. /* Return the largest (smallest) number representable in this format */
  16. fixpoint fp_max() {  return(~OVERFLOWMASK); }
  17. fixpoint fp_min() {  return(OVERFLOWMASK); }
  18.  
  19. /* return integer part */
  20. int fp_integer(fixpoint x)
  21. {
  22.   fixpoint floor = x >> LOBITS;
  23.   fixpoint bits = (x & LOMASK);
  24.  
  25.   if (x < 0) {
  26.      if (!bits) return(floor);
  27.      return(floor+1);
  28.   }
  29.   else return (floor);
  30. }
  31.  
  32. /* return fraction part, in bits */
  33. int fp_fraction(fixpoint x)
  34. {
  35.   if (x < 0) return(TwosComplement(x) & LOMASK);
  36.   else return(x & LOMASK);
  37. }
  38.  
  39. /* return fraction part, as a double */
  40. double fp_fraction_double(fixpoint x)
  41. {
  42.   if (x < 0) return(-(TwosComplement(x) & LOMASK) / ((double) (1 << LOBITS)));
  43.   else return((x & LOMASK) / ((double) (1 << LOBITS)));
  44. }
  45.  
  46. /* in 2s complement, integers have all lower order bits = 0,
  47.     and truncating the fraction = floor.
  48.  */
  49. fixpoint fp_floor(fixpoint x)
  50. {
  51.   fixpoint answer = x & HIMASK;
  52.   return (answer);
  53. }
  54.  
  55.  
  56. /* The std integer truncating divide does a floor operation,
  57.   but for negative numbers it does floor(abs(x/y)), so
  58.   we need to case on sign of x to fix that.
  59.  
  60.  if x%y has no remainder, then the divide is on an integer
  61.  
  62. Overflow can't happen in fp_floor_div, since x/y <= x for y>=1, which
  63. must be the case since y is an integer and y>0.
  64.  
  65. */
  66.  
  67. fixpoint fp_floor_div(fixpoint x, fixpoint y)
  68. {
  69.   fixpoint answer;
  70.   
  71.   if (y == 0) {
  72.      printf("Error: fp_floor_div -- can't divide by 0!\n");
  73.      return(fp_fix(0));
  74.   }
  75.   else if (y < 0) {
  76.      printf("Sorry, fp_floor_div(a,b) doesn't currently work for b<0.\n");
  77.      return(fp_fix(0));
  78.   }
  79.   else if (x >= 0) {
  80.      fixpoint tmp = x/y;
  81.      answer = tmp << LOBITS;
  82.   }
  83.   else {
  84.      fixpoint tmp = ((x/y) + (((x % y) == 0) ? 0 : -1));
  85.      answer = tmp << LOBITS;
  86.   }
  87.   return (answer);
  88. }
  89.  
  90.  
  91. fixpoint fp_multiply(fixpoint x, fixpoint y)
  92. {
  93.   fixpoint answer;
  94.   unsigned int xhi, xlo, yhi, ylo, tmp;
  95.  
  96.   fp_error = 0;
  97.  
  98.   xhi = (x<0) ? -fp_integer(x) : fp_integer(x);
  99.   yhi = (y<0) ? -fp_integer(y) : fp_integer(y);
  100.   xlo = fp_fraction(x);
  101.   ylo = fp_fraction(y);
  102.  
  103.   /* If xhi and yhi are both < 16 bits, this can't have machine overflow
  104.       (overflow of the 32bit int). But it CAN get larger than HIBITS, or
  105.       crush the sign bit, causing overflow of our fixpoint representation.
  106.     */
  107.   tmp = (xhi * yhi);
  108.   if (tmp & (((int) OVERFLOWMASK) >> LOBITS)) {
  109.      if (fp_print_error)
  110.         printf ("ERROR: fp_multiply() xhi*yhi = 0x%08x overflows by 0x%08x\n",
  111.                   tmp,
  112.                   (tmp & (((int) OVERFLOWMASK) >> LOBITS)));
  113.      fp_error = -1;
  114.      abort();
  115.      return(0);
  116.   }
  117.  
  118.   answer = (tmp << LOBITS) +
  119.            (xhi * ylo) + (xlo * yhi) +
  120.               (((xlo * ylo) >> LOBITS) & LOMASK);
  121.  
  122.   if (answer & ~(HIMASK | LOMASK)) { 
  123.      if (fp_print_error)
  124.         printf ("ERROR: fp_multiply() answer = 0x%08x, overflow 0x%08x\n",
  125.                   answer, (answer & ~(HIMASK | LOMASK)));
  126.      fp_error = -2;
  127.      abort();
  128.      return(0);
  129.   }
  130.  
  131.   /* Also, it can fill up the high order (sign) bit, which is overflow. */
  132.   if (answer & SIGNBIT) {
  133.      if (fp_print_error)
  134.         printf ("ERROR: fp_multiply() SIGNBIT overflow for answer = 0x%08x\n",
  135.                   answer);
  136.      fp_error = -3;
  137.      return(0);
  138.   }
  139.  
  140.   /* "this should never happen" */
  141.   if ((xlo >> 16) && (ylo >> 16)) {
  142.      if (fp_print_error)
  143.         printf("++++++++++++  ERROR -- OVERFLOW OF LOW BITS**********\n");
  144.      fp_error = -4;
  145.      abort();
  146.      return(0);
  147.   }
  148.  
  149.   if (((x<0) && (y>0)) || ((x>0) && (y<0))) answer = -answer;
  150.  
  151.   return(answer);
  152. }
  153.  
  154.  
  155. /* Turn a double into a fixpoint number int two's complement. 
  156.     Negative numbers become:  ~a + 1
  157.  */
  158. fixpoint fp_fix(double x)
  159. {
  160.   int negative = (x < 0);
  161.   double i;                /* integer part */
  162.   double fraction;    /* fraction part, positive only */
  163.   fixpoint p;            /* fixpoint version of abs(x) */
  164.  
  165.  
  166.   if ((x < fp_integer(fp_min())) || (x >= fp_integer(fp_max())))
  167.      printf("sorry, %g doesn't fit in (%d hi, %d lo) bit fix point.\n",
  168.               x, HIBITS, LOBITS);
  169.  
  170.   x = fabs(x);
  171.   i = floor(x);
  172.   fraction = x - i;
  173.  
  174.   p = (((int) i) << LOBITS) | ((int) ((x-i)*(1<<LOBITS)));
  175.  
  176.   if (negative) return TwosComplement(p);
  177.   else return (p);
  178. }
  179.  
  180. /* Print as a binary number */
  181. void fp_printb(fixpoint x)
  182. {
  183.   int i;
  184.   unsigned int mask;
  185.  
  186.   if (x<0) x = -x;
  187.   mask = 1 << (HIBITS + LOBITS - 1); /* start printing hi bit */
  188.   for(i=0; i<HIBITS; i++) {
  189.      if (mask & x) putchar('1');
  190.      else          putchar('0');
  191.      mask = mask >> 1;
  192.   }
  193.   putchar('.');
  194.   for(i=0; i<LOBITS; i++) {
  195.      if (mask & x) putchar('1');
  196.      else          putchar('0');
  197.      mask = mask >> 1;
  198.   }
  199. }
  200.  
  201. /* Print as a hex number, but gets things a bit screwed
  202.     up unless HIBITS=LOBITS=16
  203. */
  204. void fp_printx(fixpoint x)
  205. {
  206.   printf("%4x.%04x", fp_integer(x) & LOMASK, fp_fraction(x));
  207. }
  208.  
  209. double fp_double(fixpoint x)
  210. {
  211.   return(((double) fp_integer(x)) + fp_fraction_double(x));
  212. }
  213.  
  214. void fp_print(fixpoint x)
  215. {
  216.   printf("%.11g", fp_double(x));
  217. }
  218.  
  219.  
  220. void printnbits()
  221. {
  222.   printf("HIBITS = %d, LOBITS = %d, OVERFLOWMASK = 0x%8x, SIGNBIT = 0x%08x\n",
  223.             HIBITS, LOBITS, OVERFLOWMASK, SIGNBIT);
  224.   printf("HIMASK = 0x%08x, LOMASK = 0x%08x\n", HIMASK, LOMASK);
  225.   printf("overflowmask >> lobits = 0x%08x\n", ((int) OVERFLOWMASK) >> LOBITS);
  226. }
  227.